home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / blas / blas.c
Encoding:
C/C++ Source or Header  |  2002-04-18  |  54.1 KB  |  2,183 lines
  1. /* blas/blas.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 Gerard Jungman & Brian 
  4.  * Gough
  5.  * 
  6.  * This program is free software; you can redistribute it and/or modify
  7.  * it under the terms of the GNU General Public License as published by
  8.  * the Free Software Foundation; either version 2 of the License, or (at
  9.  * your option) any later version.
  10.  * 
  11.  * This program is distributed in the hope that it will be useful, but
  12.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14.  * General Public License for more details.
  15.  * 
  16.  * You should have received a copy of the GNU General Public License
  17.  * along with this program; if not, write to the Free Software
  18.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  */
  20.  
  21. /* GSL implementation of BLAS operations for vectors and dense
  22.  * matrices.  Note that GSL native storage is row-major.  */
  23.  
  24. #include <config.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_errno.h>
  27. #include <gsl/gsl_cblas.h>
  28. #include <gsl/gsl_cblas.h>
  29. #include <gsl/gsl_blas_types.h>
  30. #include <gsl/gsl_blas.h>
  31.  
  32. /* ========================================================================
  33.  * Level 1
  34.  * ========================================================================
  35.  */
  36.  
  37. /* CBLAS defines vector sizes in terms of int. GSL defines sizes in
  38.    terms of size_t, so we need to convert these into integers.  There
  39.    is the possibility of overflow here. FIXME: Maybe this could be
  40.    caught */
  41.  
  42. #define INT(X) ((int)(X))
  43.  
  44. int
  45. gsl_blas_sdsdot (float alpha, const gsl_vector_float * X,
  46.          const gsl_vector_float * Y, float *result)
  47. {
  48.   if (X->size == Y->size)
  49.     {
  50.       *result =
  51.     cblas_sdsdot (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
  52.               INT (Y->stride));
  53.       return GSL_SUCCESS;
  54.     }
  55.   else
  56.     {
  57.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  58.     }
  59. }
  60.  
  61. int
  62. gsl_blas_dsdot (const gsl_vector_float * X, const gsl_vector_float * Y,
  63.         double *result)
  64. {
  65.   if (X->size == Y->size)
  66.     {
  67.       *result =
  68.     cblas_dsdot (INT (X->size), X->data, INT (X->stride), Y->data,
  69.              INT (Y->stride));
  70.       return GSL_SUCCESS;
  71.     }
  72.   else
  73.     {
  74.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  75.     }
  76. }
  77.  
  78. int
  79. gsl_blas_sdot (const gsl_vector_float * X, const gsl_vector_float * Y,
  80.            float *result)
  81. {
  82.   if (X->size == Y->size)
  83.     {
  84.       *result =
  85.     cblas_sdot (INT (X->size), X->data, INT (X->stride), Y->data,
  86.             INT (Y->stride));
  87.       return GSL_SUCCESS;
  88.     }
  89.   else
  90.     {
  91.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  92.     }
  93. }
  94.  
  95. int
  96. gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result)
  97. {
  98.   if (X->size == Y->size)
  99.     {
  100.       *result =
  101.     cblas_ddot (INT (X->size), X->data, INT (X->stride), Y->data,
  102.             INT (Y->stride));
  103.       return GSL_SUCCESS;
  104.     }
  105.   else
  106.     {
  107.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  108.     }
  109. }
  110.  
  111.  
  112. int
  113. gsl_blas_cdotu (const gsl_vector_complex_float * X,
  114.         const gsl_vector_complex_float * Y, gsl_complex_float * dotu)
  115. {
  116.   if (X->size == Y->size)
  117.     {
  118.       cblas_cdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
  119.                INT (Y->stride), GSL_COMPLEX_P (dotu));
  120.       return GSL_SUCCESS;
  121.     }
  122.   else
  123.     {
  124.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  125.     }
  126. }
  127.  
  128.  
  129. int
  130. gsl_blas_cdotc (const gsl_vector_complex_float * X,
  131.         const gsl_vector_complex_float * Y, gsl_complex_float * dotc)
  132. {
  133.   if (X->size == Y->size)
  134.     {
  135.       cblas_cdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
  136.                INT (Y->stride), GSL_COMPLEX_P (dotc));
  137.       return GSL_SUCCESS;
  138.     }
  139.   else
  140.     {
  141.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  142.     }
  143. }
  144.  
  145.  
  146. int
  147. gsl_blas_zdotu (const gsl_vector_complex * X, const gsl_vector_complex * Y,
  148.         gsl_complex * dotu)
  149. {
  150.   if (X->size == Y->size)
  151.     {
  152.       cblas_zdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
  153.                INT (Y->stride), GSL_COMPLEX_P (dotu));
  154.       return GSL_SUCCESS;
  155.     }
  156.   else
  157.     {
  158.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  159.     }
  160. }
  161.  
  162.  
  163. int
  164. gsl_blas_zdotc (const gsl_vector_complex * X, const gsl_vector_complex * Y,
  165.         gsl_complex * dotc)
  166. {
  167.   if (X->size == Y->size)
  168.     {
  169.       cblas_zdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
  170.                INT (Y->stride), GSL_COMPLEX_P (dotc));
  171.       return GSL_SUCCESS;
  172.     }
  173.   else
  174.     {
  175.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  176.     }
  177. }
  178.  
  179. /* Norms of vectors */
  180.  
  181. float
  182. gsl_blas_snrm2 (const gsl_vector_float * X)
  183. {
  184.   return cblas_snrm2 (INT (X->size), X->data, INT (X->stride));
  185. }
  186.  
  187. double
  188. gsl_blas_dnrm2 (const gsl_vector * X)
  189. {
  190.   return cblas_dnrm2 (INT (X->size), X->data, INT (X->stride));
  191. }
  192.  
  193. float
  194. gsl_blas_scnrm2 (const gsl_vector_complex_float * X)
  195. {
  196.   return cblas_scnrm2 (INT (X->size), X->data, INT (X->stride));
  197. }
  198.  
  199. double
  200. gsl_blas_dznrm2 (const gsl_vector_complex * X)
  201. {
  202.   return cblas_dznrm2 (INT (X->size), X->data, INT (X->stride));
  203. }
  204.  
  205. /* Absolute sums of vectors */
  206.  
  207. float
  208. gsl_blas_sasum (const gsl_vector_float * X)
  209. {
  210.   return cblas_sasum (INT (X->size), X->data, INT (X->stride));
  211. }
  212.  
  213. double
  214. gsl_blas_dasum (const gsl_vector * X)
  215. {
  216.   return cblas_dasum (INT (X->size), X->data, INT (X->stride));
  217. }
  218.  
  219. float
  220. gsl_blas_scasum (const gsl_vector_complex_float * X)
  221. {
  222.   return cblas_scasum (INT (X->size), X->data, INT (X->stride));
  223. }
  224.  
  225. double
  226. gsl_blas_dzasum (const gsl_vector_complex * X)
  227. {
  228.   return cblas_dzasum (INT (X->size), X->data, INT (X->stride));
  229. }
  230.  
  231. /* Maximum elements of vectors */
  232.  
  233. CBLAS_INDEX_t
  234. gsl_blas_isamax (const gsl_vector_float * X)
  235. {
  236.   return cblas_isamax (INT (X->size), X->data, INT (X->stride));
  237. }
  238.  
  239. CBLAS_INDEX_t
  240. gsl_blas_idamax (const gsl_vector * X)
  241. {
  242.   return cblas_idamax (INT (X->size), X->data, INT (X->stride));
  243. }
  244.  
  245. CBLAS_INDEX_t
  246. gsl_blas_icamax (const gsl_vector_complex_float * X)
  247. {
  248.   return cblas_icamax (INT (X->size), X->data, INT (X->stride));
  249. }
  250.  
  251. CBLAS_INDEX_t
  252. gsl_blas_izamax (const gsl_vector_complex * X)
  253. {
  254.   return cblas_izamax (INT (X->size), X->data, INT (X->stride));
  255. }
  256.  
  257.  
  258. /* Swap vectors */
  259.  
  260. int
  261. gsl_blas_sswap (gsl_vector_float * X, gsl_vector_float * Y)
  262. {
  263.   if (X->size == Y->size)
  264.     {
  265.       cblas_sswap (INT (X->size), X->data, INT (X->stride), Y->data,
  266.            INT (Y->stride));
  267.       return GSL_SUCCESS;
  268.     }
  269.   else
  270.     {
  271.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  272.     }
  273. }
  274.  
  275. int
  276. gsl_blas_dswap (gsl_vector * X, gsl_vector * Y)
  277. {
  278.   if (X->size == Y->size)
  279.     {
  280.       cblas_dswap (INT (X->size), X->data, INT (X->stride), Y->data,
  281.            INT (Y->stride));
  282.       return GSL_SUCCESS;
  283.     }
  284.   else
  285.     {
  286.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  287.     };
  288. }
  289.  
  290. int
  291. gsl_blas_cswap (gsl_vector_complex_float * X, gsl_vector_complex_float * Y)
  292. {
  293.   if (X->size == Y->size)
  294.     {
  295.       cblas_cswap (INT (X->size), X->data, INT (X->stride), Y->data,
  296.            INT (Y->stride));
  297.       return GSL_SUCCESS;
  298.     }
  299.   else
  300.     {
  301.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  302.     }
  303. }
  304.  
  305. int
  306. gsl_blas_zswap (gsl_vector_complex * X, gsl_vector_complex * Y)
  307. {
  308.   if (X->size == Y->size)
  309.     {
  310.       cblas_zswap (INT (X->size), X->data, INT (X->stride), Y->data,
  311.            INT (Y->stride));
  312.       return GSL_SUCCESS;
  313.     }
  314.   else
  315.     {
  316.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  317.     }
  318. }
  319.  
  320.  
  321. /* Copy vectors */
  322.  
  323. int
  324. gsl_blas_scopy (const gsl_vector_float * X, gsl_vector_float * Y)
  325. {
  326.   if (X->size == Y->size)
  327.     {
  328.       cblas_scopy (INT (X->size), X->data, INT (X->stride), Y->data,
  329.            INT (Y->stride));
  330.       return GSL_SUCCESS;
  331.     }
  332.   else
  333.     {
  334.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  335.     }
  336. }
  337.  
  338. int
  339. gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y)
  340. {
  341.   if (X->size == Y->size)
  342.     {
  343.       cblas_dcopy (INT (X->size), X->data, INT (X->stride), Y->data,
  344.            INT (Y->stride));
  345.       return GSL_SUCCESS;
  346.     }
  347.   else
  348.     {
  349.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  350.     }
  351. }
  352.  
  353. int
  354. gsl_blas_ccopy (const gsl_vector_complex_float * X,
  355.         gsl_vector_complex_float * Y)
  356. {
  357.   if (X->size == Y->size)
  358.     {
  359.       cblas_ccopy (INT (X->size), X->data, INT (X->stride), Y->data,
  360.            INT (Y->stride));
  361.       return GSL_SUCCESS;
  362.     }
  363.   else
  364.     {
  365.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  366.     }
  367. }
  368.  
  369. int
  370. gsl_blas_zcopy (const gsl_vector_complex * X, gsl_vector_complex * Y)
  371. {
  372.   if (X->size == Y->size)
  373.     {
  374.       cblas_zcopy (INT (X->size), X->data, INT (X->stride), Y->data,
  375.            INT (Y->stride));
  376.       return GSL_SUCCESS;
  377.     }
  378.   else
  379.     {
  380.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  381.     }
  382. }
  383.  
  384.  
  385. /* Compute Y = alpha X + Y */
  386.  
  387. int
  388. gsl_blas_saxpy (float alpha, const gsl_vector_float * X, gsl_vector_float * Y)
  389. {
  390.   if (X->size == Y->size)
  391.     {
  392.       cblas_saxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
  393.            INT (Y->stride));
  394.       return GSL_SUCCESS;
  395.     }
  396.   else
  397.     {
  398.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  399.     }
  400. }
  401.  
  402. int
  403. gsl_blas_daxpy (double alpha, const gsl_vector * X, gsl_vector * Y)
  404. {
  405.   if (X->size == Y->size)
  406.     {
  407.       cblas_daxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
  408.            INT (Y->stride));
  409.       return GSL_SUCCESS;
  410.     }
  411.   else
  412.     {
  413.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  414.     }
  415. }
  416.  
  417. int
  418. gsl_blas_caxpy (const gsl_complex_float alpha,
  419.         const gsl_vector_complex_float * X,
  420.         gsl_vector_complex_float * Y)
  421. {
  422.   if (X->size == Y->size)
  423.     {
  424.       cblas_caxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
  425.            INT (X->stride), Y->data, INT (Y->stride));
  426.       return GSL_SUCCESS;
  427.     }
  428.   else
  429.     {
  430.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  431.     }
  432. }
  433.  
  434. int
  435. gsl_blas_zaxpy (const gsl_complex alpha, const gsl_vector_complex * X,
  436.         gsl_vector_complex * Y)
  437. {
  438.   if (X->size == Y->size)
  439.     {
  440.       cblas_zaxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
  441.            INT (X->stride), Y->data, INT (Y->stride));
  442.       return GSL_SUCCESS;
  443.     }
  444.   else
  445.     {
  446.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  447.     }
  448. }
  449.  
  450. /* Generate rotation */
  451.  
  452. int
  453. gsl_blas_srotg (float a[], float b[], float c[], float s[])
  454. {
  455.   cblas_srotg (a, b, c, s);
  456.   return GSL_SUCCESS;
  457. }
  458.  
  459. int
  460. gsl_blas_drotg (double a[], double b[], double c[], double s[])
  461. {
  462.   cblas_drotg (a, b, c, s);
  463.   return GSL_SUCCESS;
  464. }
  465.  
  466. /* Apply rotation to vectors */
  467.  
  468. int
  469. gsl_blas_srot (gsl_vector_float * X, gsl_vector_float * Y, float c, float s)
  470. {
  471.   if (X->size == Y->size)
  472.     {
  473.       cblas_srot (INT (X->size), X->data, INT (X->stride), Y->data,
  474.           INT (Y->stride), c, s);
  475.       return GSL_SUCCESS;
  476.     }
  477.   else
  478.     {
  479.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  480.     }
  481. }
  482.  
  483. int
  484. gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const double c, const double s)
  485. {
  486.   if (X->size == Y->size)
  487.     {
  488.       cblas_drot (INT (X->size), X->data, INT (X->stride), Y->data,
  489.           INT (Y->stride), c, s);
  490.       return GSL_SUCCESS;
  491.     }
  492.   else
  493.     {
  494.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  495.     }
  496. }
  497.  
  498.  
  499. /* Generate modified rotation */
  500.  
  501. int
  502. gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[])
  503. {
  504.   cblas_srotmg (d1, d2, b1, b2, P);
  505.   return GSL_SUCCESS;
  506. }
  507.  
  508. int
  509. gsl_blas_drotmg (double d1[], double d2[], double b1[], double b2, double P[])
  510. {
  511.   cblas_drotmg (d1, d2, b1, b2, P);
  512.   return GSL_SUCCESS;
  513. }
  514.  
  515.  
  516. /* Apply modified rotation */
  517.  
  518. int
  519. gsl_blas_srotm (gsl_vector_float * X, gsl_vector_float * Y, const float P[])
  520. {
  521.   if (X->size == Y->size)
  522.     {
  523.       cblas_srotm (INT (X->size), X->data, INT (X->stride), Y->data,
  524.            INT (Y->stride), P);
  525.       return GSL_SUCCESS;
  526.     }
  527.   else
  528.     {
  529.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  530.     }
  531. }
  532.  
  533. int
  534. gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[])
  535. {
  536.   if (X->size != Y->size)
  537.     {
  538.       cblas_drotm (INT (X->size), X->data, INT (X->stride), Y->data,
  539.            INT (Y->stride), P);
  540.       return GSL_SUCCESS;
  541.     }
  542.   else
  543.     {
  544.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  545.     }
  546. }
  547.  
  548.  
  549. /* Scale vector */
  550.  
  551. void
  552. gsl_blas_sscal (float alpha, gsl_vector_float * X)
  553. {
  554.   cblas_sscal (INT (X->size), alpha, X->data, INT (X->stride));
  555. }
  556.  
  557. void
  558. gsl_blas_dscal (double alpha, gsl_vector * X)
  559. {
  560.   cblas_dscal (INT (X->size), alpha, X->data, INT (X->stride));
  561. }
  562.  
  563. void
  564. gsl_blas_cscal (const gsl_complex_float alpha, gsl_vector_complex_float * X)
  565. {
  566.   cblas_cscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
  567.            INT (X->stride));
  568. }
  569.  
  570. void
  571. gsl_blas_zscal (const gsl_complex alpha, gsl_vector_complex * X)
  572. {
  573.   cblas_zscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
  574.            INT (X->stride));
  575. }
  576.  
  577. void
  578. gsl_blas_csscal (float alpha, gsl_vector_complex_float * X)
  579. {
  580.   cblas_csscal (INT (X->size), alpha, X->data, INT (X->stride));
  581. }
  582.  
  583. void
  584. gsl_blas_zdscal (double alpha, gsl_vector_complex * X)
  585. {
  586.   cblas_zdscal (INT (X->size), alpha, X->data, INT (X->stride));
  587. }
  588.  
  589. /* ===========================================================================
  590.  * Level 2
  591.  * ===========================================================================
  592.  */
  593.  
  594. /* GEMV */
  595.  
  596. int
  597. gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA, float alpha,
  598.         const gsl_matrix_float * A, const gsl_vector_float * X,
  599.         float beta, gsl_vector_float * Y)
  600. {
  601.   const size_t M = A->size1;
  602.   const size_t N = A->size2;
  603.  
  604.   if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
  605.       || (TransA == CblasTrans && M == X->size && N == Y->size))
  606.     {
  607.       cblas_sgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
  608.            INT (A->tda), X->data, INT (X->stride), beta, Y->data,
  609.            INT (Y->stride));
  610.       return GSL_SUCCESS;
  611.     }
  612.   else
  613.     {
  614.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  615.     }
  616. }
  617.  
  618.  
  619. int
  620. gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,
  621.         const gsl_vector * X, double beta, gsl_vector * Y)
  622. {
  623.   const size_t M = A->size1;
  624.   const size_t N = A->size2;
  625.  
  626.   if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
  627.       || (TransA == CblasTrans && M == X->size && N == Y->size))
  628.     {
  629.       cblas_dgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
  630.            INT (A->tda), X->data, INT (X->stride), beta, Y->data,
  631.            INT (Y->stride));
  632.       return GSL_SUCCESS;
  633.     }
  634.   else
  635.     {
  636.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  637.     }
  638. }
  639.  
  640.  
  641. int
  642. gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex_float alpha,
  643.         const gsl_matrix_complex_float * A,
  644.         const gsl_vector_complex_float * X,
  645.         const gsl_complex_float beta, gsl_vector_complex_float * Y)
  646. {
  647.   const size_t M = A->size1;
  648.   const size_t N = A->size2;
  649.  
  650.   if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
  651.       || (TransA == CblasTrans && M == X->size && N == Y->size))
  652.     {
  653.       cblas_cgemv (CblasRowMajor, TransA, INT (M), INT (N),
  654.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
  655.            INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
  656.            INT (Y->stride));
  657.       return GSL_SUCCESS;
  658.     }
  659.   else
  660.     {
  661.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  662.     }
  663. }
  664.  
  665.  
  666. int
  667. gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex alpha,
  668.         const gsl_matrix_complex * A, const gsl_vector_complex * X,
  669.         const gsl_complex beta, gsl_vector_complex * Y)
  670. {
  671.   const size_t M = A->size1;
  672.   const size_t N = A->size2;
  673.  
  674.   if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
  675.       || (TransA == CblasTrans && M == X->size && N == Y->size))
  676.     {
  677.       cblas_zgemv (CblasRowMajor, TransA, INT (M), INT (N),
  678.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
  679.            INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
  680.            INT (Y->stride));
  681.       return GSL_SUCCESS;
  682.     }
  683.   else
  684.     {
  685.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  686.     }
  687. }
  688.  
  689.  
  690.  
  691. /* HEMV */
  692.  
  693. int
  694. gsl_blas_chemv (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
  695.         const gsl_matrix_complex_float * A,
  696.         const gsl_vector_complex_float * X,
  697.         const gsl_complex_float beta, gsl_vector_complex_float * Y)
  698. {
  699.   const size_t M = A->size1;
  700.   const size_t N = A->size2;
  701.  
  702.   if (M != N)
  703.     {
  704.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  705.     }
  706.   else if (N != X->size || N != Y->size)
  707.     {
  708.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  709.     }
  710.  
  711.   cblas_chemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
  712.            INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
  713.            Y->data, INT (Y->stride));
  714.   return GSL_SUCCESS;
  715. }
  716.  
  717. int
  718. gsl_blas_zhemv (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
  719.         const gsl_matrix_complex * A, const gsl_vector_complex * X,
  720.         const gsl_complex beta, gsl_vector_complex * Y)
  721. {
  722.   const size_t M = A->size1;
  723.   const size_t N = A->size2;
  724.  
  725.   if (M != N)
  726.     {
  727.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  728.     }
  729.   else if (N != X->size || N != Y->size)
  730.     {
  731.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  732.     }
  733.  
  734.   cblas_zhemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
  735.            INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
  736.            Y->data, INT (Y->stride));
  737.   return GSL_SUCCESS;
  738. }
  739.  
  740.  
  741. /* SYMV */
  742.  
  743. int
  744. gsl_blas_ssymv (CBLAS_UPLO_t Uplo, float alpha, const gsl_matrix_float * A,
  745.         const gsl_vector_float * X, float beta, gsl_vector_float * Y)
  746. {
  747.   const size_t M = A->size1;
  748.   const size_t N = A->size2;
  749.  
  750.   if (M != N)
  751.     {
  752.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  753.     }
  754.   else if (N != X->size || N != Y->size)
  755.     {
  756.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  757.     }
  758.  
  759.   cblas_ssymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
  760.            X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
  761.   return GSL_SUCCESS;
  762. }
  763.  
  764. int
  765. gsl_blas_dsymv (CBLAS_UPLO_t Uplo, double alpha, const gsl_matrix * A,
  766.         const gsl_vector * X, double beta, gsl_vector * Y)
  767. {
  768.   const size_t M = A->size1;
  769.   const size_t N = A->size2;
  770.  
  771.   if (M != N)
  772.     {
  773.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  774.     }
  775.   else if (N != X->size || N != Y->size)
  776.     {
  777.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  778.     }
  779.  
  780.   cblas_dsymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
  781.            X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
  782.   return GSL_SUCCESS;
  783. }
  784.  
  785.  
  786. /* TRMV */
  787.  
  788. int
  789. gsl_blas_strmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  790.         CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
  791.         gsl_vector_float * X)
  792. {
  793.   const size_t M = A->size1;
  794.   const size_t N = A->size2;
  795.  
  796.   if (M != N)
  797.     {
  798.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  799.     }
  800.   else if (N != X->size)
  801.     {
  802.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  803.     }
  804.  
  805.   cblas_strmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  806.            INT (A->tda), X->data, INT (X->stride));
  807.   return GSL_SUCCESS;
  808. }
  809.  
  810.  
  811. int
  812. gsl_blas_dtrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  813.         CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
  814. {
  815.   const size_t M = A->size1;
  816.   const size_t N = A->size2;
  817.  
  818.   if (M != N)
  819.     {
  820.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  821.     }
  822.   else if (N != X->size)
  823.     {
  824.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  825.     }
  826.  
  827.   cblas_dtrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  828.            INT (A->tda), X->data, INT (X->stride));
  829.   return GSL_SUCCESS;
  830. }
  831.  
  832.  
  833. int
  834. gsl_blas_ctrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  835.         CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
  836.         gsl_vector_complex_float * X)
  837. {
  838.   const size_t M = A->size1;
  839.   const size_t N = A->size2;
  840.  
  841.   if (M != N)
  842.     {
  843.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  844.     }
  845.   else if (N != X->size)
  846.     {
  847.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  848.     }
  849.  
  850.   cblas_ctrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  851.            INT (A->tda), X->data, INT (X->stride));
  852.   return GSL_SUCCESS;
  853. }
  854.  
  855.  
  856. int
  857. gsl_blas_ztrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  858.         CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
  859.         gsl_vector_complex * X)
  860. {
  861.   const size_t M = A->size1;
  862.   const size_t N = A->size2;
  863.  
  864.   if (M != N)
  865.     {
  866.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  867.     }
  868.   else if (N != X->size)
  869.     {
  870.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  871.     }
  872.  
  873.   cblas_ztrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  874.            INT (A->tda), X->data, INT (X->stride));
  875.   return GSL_SUCCESS;
  876. }
  877.  
  878.  
  879. /* TRSV */
  880.  
  881. int
  882. gsl_blas_strsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  883.         CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
  884.         gsl_vector_float * X)
  885. {
  886.   const size_t M = A->size1;
  887.   const size_t N = A->size2;
  888.  
  889.   if (M != N)
  890.     {
  891.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  892.     }
  893.   else if (N != X->size)
  894.     {
  895.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  896.     }
  897.  
  898.   cblas_strsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  899.            INT (A->tda), X->data, INT (X->stride));
  900.   return GSL_SUCCESS;
  901. }
  902.  
  903.  
  904. int
  905. gsl_blas_dtrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  906.         CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
  907. {
  908.   const size_t M = A->size1;
  909.   const size_t N = A->size2;
  910.  
  911.   if (M != N)
  912.     {
  913.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  914.     }
  915.   else if (N != X->size)
  916.     {
  917.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  918.     }
  919.  
  920.   cblas_dtrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  921.            INT (A->tda), X->data, INT (X->stride));
  922.   return GSL_SUCCESS;
  923. }
  924.  
  925.  
  926. int
  927. gsl_blas_ctrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  928.         CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
  929.         gsl_vector_complex_float * X)
  930. {
  931.   const size_t M = A->size1;
  932.   const size_t N = A->size2;
  933.  
  934.   if (M != N)
  935.     {
  936.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  937.     }
  938.   else if (N != X->size)
  939.     {
  940.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  941.     }
  942.  
  943.   cblas_ctrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  944.            INT (A->tda), X->data, INT (X->stride));
  945.   return GSL_SUCCESS;
  946. }
  947.  
  948.  
  949. int
  950. gsl_blas_ztrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
  951.         CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
  952.         gsl_vector_complex * X)
  953. {
  954.   const size_t M = A->size1;
  955.   const size_t N = A->size2;
  956.  
  957.   if (M != N)
  958.     {
  959.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  960.     }
  961.   else if (N != X->size)
  962.     {
  963.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  964.     }
  965.  
  966.   cblas_ztrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
  967.            INT (A->tda), X->data, INT (X->stride));
  968.   return GSL_SUCCESS;
  969. }
  970.  
  971.  
  972. /* GER */
  973.  
  974. int
  975. gsl_blas_sger (float alpha, const gsl_vector_float * X,
  976.            const gsl_vector_float * Y, gsl_matrix_float * A)
  977. {
  978.   const size_t M = A->size1;
  979.   const size_t N = A->size2;
  980.  
  981.   if (X->size == M && Y->size == N)
  982.     {
  983.       cblas_sger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
  984.           INT (X->stride), Y->data, INT (Y->stride), A->data,
  985.           INT (A->tda));
  986.       return GSL_SUCCESS;
  987.     }
  988.   else
  989.     {
  990.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  991.     }
  992. }
  993.  
  994.  
  995. int
  996. gsl_blas_dger (double alpha, const gsl_vector * X, const gsl_vector * Y,
  997.            gsl_matrix * A)
  998. {
  999.   const size_t M = A->size1;
  1000.   const size_t N = A->size2;
  1001.  
  1002.   if (X->size == M && Y->size == N)
  1003.     {
  1004.       cblas_dger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
  1005.           INT (X->stride), Y->data, INT (Y->stride), A->data,
  1006.           INT (A->tda));
  1007.       return GSL_SUCCESS;
  1008.     }
  1009.   else
  1010.     {
  1011.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1012.     }
  1013. }
  1014.  
  1015.  
  1016. /* GERU */
  1017.  
  1018. int
  1019. gsl_blas_cgeru (const gsl_complex_float alpha,
  1020.         const gsl_vector_complex_float * X,
  1021.         const gsl_vector_complex_float * Y,
  1022.         gsl_matrix_complex_float * A)
  1023. {
  1024.   const size_t M = A->size1;
  1025.   const size_t N = A->size2;
  1026.  
  1027.   if (X->size == M && Y->size == N)
  1028.     {
  1029.       cblas_cgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
  1030.            X->data, INT (X->stride), Y->data, INT (Y->stride),
  1031.            A->data, INT (A->tda));
  1032.       return GSL_SUCCESS;
  1033.     }
  1034.   else
  1035.     {
  1036.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1037.     }
  1038. }
  1039.  
  1040. int
  1041. gsl_blas_zgeru (const gsl_complex alpha, const gsl_vector_complex * X,
  1042.         const gsl_vector_complex * Y, gsl_matrix_complex * A)
  1043. {
  1044.   const size_t M = A->size1;
  1045.   const size_t N = A->size2;
  1046.  
  1047.   if (X->size == M && Y->size == N)
  1048.     {
  1049.       cblas_zgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
  1050.            X->data, INT (X->stride), Y->data, INT (Y->stride),
  1051.            A->data, INT (A->tda));
  1052.       return GSL_SUCCESS;
  1053.     }
  1054.   else
  1055.     {
  1056.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1057.     }
  1058. }
  1059.  
  1060.  
  1061. /* GERC */
  1062.  
  1063. int
  1064. gsl_blas_cgerc (const gsl_complex_float alpha,
  1065.         const gsl_vector_complex_float * X,
  1066.         const gsl_vector_complex_float * Y,
  1067.         gsl_matrix_complex_float * A)
  1068. {
  1069.   const size_t M = A->size1;
  1070.   const size_t N = A->size2;
  1071.  
  1072.   if (X->size == M && Y->size == N)
  1073.     {
  1074.       cblas_cgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
  1075.            X->data, INT (X->stride), Y->data, INT (Y->stride),
  1076.            A->data, INT (A->tda));
  1077.       return GSL_SUCCESS;
  1078.     }
  1079.   else
  1080.     {
  1081.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1082.     }
  1083. }
  1084.  
  1085.  
  1086. int
  1087. gsl_blas_zgerc (const gsl_complex alpha, const gsl_vector_complex * X,
  1088.         const gsl_vector_complex * Y, gsl_matrix_complex * A)
  1089. {
  1090.   const size_t M = A->size1;
  1091.   const size_t N = A->size2;
  1092.  
  1093.   if (X->size == M && Y->size == N)
  1094.     {
  1095.       cblas_zgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
  1096.            X->data, INT (X->stride), Y->data, INT (Y->stride),
  1097.            A->data, INT (A->tda));
  1098.       return GSL_SUCCESS;
  1099.     }
  1100.   else
  1101.     {
  1102.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1103.     }
  1104. }
  1105.  
  1106. /* HER */
  1107.  
  1108. int
  1109. gsl_blas_cher (CBLAS_UPLO_t Uplo, float alpha,
  1110.            const gsl_vector_complex_float * X,
  1111.            gsl_matrix_complex_float * A)
  1112. {
  1113.   const size_t M = A->size1;
  1114.   const size_t N = A->size2;
  1115.  
  1116.   if (M != N)
  1117.     {
  1118.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1119.     }
  1120.   else if (X->size != N)
  1121.     {
  1122.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1123.     }
  1124.  
  1125.   cblas_cher (CblasRowMajor, Uplo, INT (M), alpha, X->data, INT (X->stride),
  1126.           A->data, INT (A->tda));
  1127.   return GSL_SUCCESS;
  1128. }
  1129.  
  1130.  
  1131. int
  1132. gsl_blas_zher (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector_complex * X,
  1133.            gsl_matrix_complex * A)
  1134. {
  1135.   const size_t M = A->size1;
  1136.   const size_t N = A->size2;
  1137.  
  1138.   if (M != N)
  1139.     {
  1140.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1141.     }
  1142.   else if (X->size != N)
  1143.     {
  1144.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1145.     }
  1146.  
  1147.   cblas_zher (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
  1148.           A->data, INT (A->tda));
  1149.   return GSL_SUCCESS;
  1150. }
  1151.  
  1152.  
  1153. /* HER2 */
  1154.  
  1155. int
  1156. gsl_blas_cher2 (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
  1157.         const gsl_vector_complex_float * X,
  1158.         const gsl_vector_complex_float * Y,
  1159.         gsl_matrix_complex_float * A)
  1160. {
  1161.   const size_t M = A->size1;
  1162.   const size_t N = A->size2;
  1163.  
  1164.   if (M != N)
  1165.     {
  1166.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1167.     }
  1168.   else if (X->size != N || Y->size != N)
  1169.     {
  1170.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1171.     }
  1172.  
  1173.   cblas_cher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
  1174.            INT (X->stride), Y->data, INT (Y->stride), A->data,
  1175.            INT (A->tda));
  1176.   return GSL_SUCCESS;
  1177. }
  1178.  
  1179.  
  1180. int
  1181. gsl_blas_zher2 (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
  1182.         const gsl_vector_complex * X, const gsl_vector_complex * Y,
  1183.         gsl_matrix_complex * A)
  1184. {
  1185.   const size_t M = A->size1;
  1186.   const size_t N = A->size2;
  1187.  
  1188.   if (M != N)
  1189.     {
  1190.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1191.     }
  1192.   else if (X->size != N || Y->size != N)
  1193.     {
  1194.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1195.     }
  1196.  
  1197.   cblas_zher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
  1198.            INT (X->stride), Y->data, INT (Y->stride), A->data,
  1199.            INT (A->tda));
  1200.   return GSL_SUCCESS;
  1201. }
  1202.  
  1203.  
  1204. /* SYR */
  1205.  
  1206. int
  1207. gsl_blas_ssyr (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
  1208.            gsl_matrix_float * A)
  1209. {
  1210.   const size_t M = A->size1;
  1211.   const size_t N = A->size2;
  1212.  
  1213.   if (M != N)
  1214.     {
  1215.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1216.     }
  1217.   else if (X->size != N)
  1218.     {
  1219.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1220.     }
  1221.  
  1222.   cblas_ssyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
  1223.           A->data, INT (A->tda));
  1224.   return GSL_SUCCESS;
  1225. }
  1226.  
  1227.  
  1228. int
  1229. gsl_blas_dsyr (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
  1230.            gsl_matrix * A)
  1231. {
  1232.   const size_t M = A->size1;
  1233.   const size_t N = A->size2;
  1234.  
  1235.   if (M != N)
  1236.     {
  1237.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1238.     }
  1239.   else if (X->size != N)
  1240.     {
  1241.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1242.     }
  1243.  
  1244.   cblas_dsyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
  1245.           A->data, INT (A->tda));
  1246.   return GSL_SUCCESS;
  1247. }
  1248.  
  1249.  
  1250. /* SYR2 */
  1251.  
  1252. int
  1253. gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
  1254.         const gsl_vector_float * Y, gsl_matrix_float * A)
  1255. {
  1256.   const size_t M = A->size1;
  1257.   const size_t N = A->size2;
  1258.  
  1259.   if (M != N)
  1260.     {
  1261.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1262.     }
  1263.   else if (X->size != N || Y->size != N)
  1264.     {
  1265.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1266.     }
  1267.  
  1268.   cblas_ssyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
  1269.            Y->data, INT (Y->stride), A->data, INT (A->tda));
  1270.   return GSL_SUCCESS;
  1271. }
  1272.  
  1273.  
  1274. int
  1275. gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
  1276.         const gsl_vector * Y, gsl_matrix * A)
  1277. {
  1278.   const size_t M = A->size1;
  1279.   const size_t N = A->size2;
  1280.  
  1281.   if (M != N)
  1282.     {
  1283.       GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
  1284.     }
  1285.   else if (X->size != N || Y->size != N)
  1286.     {
  1287.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1288.     }
  1289.  
  1290.   cblas_dsyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
  1291.            Y->data, INT (Y->stride), A->data, INT (A->tda));
  1292.   return GSL_SUCCESS;
  1293. }
  1294.  
  1295.  
  1296. /*
  1297.  * ===========================================================================
  1298.  * Prototypes for level 3 BLAS
  1299.  * ===========================================================================
  1300.  */
  1301.  
  1302.  
  1303. /* GEMM */
  1304.  
  1305. int
  1306. gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
  1307.         float alpha, const gsl_matrix_float * A,
  1308.         const gsl_matrix_float * B, float beta, gsl_matrix_float * C)
  1309. {
  1310.   const size_t M = C->size1;
  1311.   const size_t N = C->size2;
  1312.   const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  1313.   const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  1314.   const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  1315.   const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
  1316.  
  1317.   if (M == MA && N == NB && NA == MB)    /* [MxN] = [MAxNA][MBxNB] */
  1318.     {
  1319.       cblas_sgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
  1320.            alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
  1321.            C->data, INT (C->tda));
  1322.       return GSL_SUCCESS;
  1323.     }
  1324.   else
  1325.     {
  1326.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1327.     }
  1328. }
  1329.  
  1330.  
  1331. int
  1332. gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
  1333.         double alpha, const gsl_matrix * A, const gsl_matrix * B,
  1334.         double beta, gsl_matrix * C)
  1335. {
  1336.   const size_t M = C->size1;
  1337.   const size_t N = C->size2;
  1338.   const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  1339.   const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  1340.   const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  1341.   const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
  1342.  
  1343.   if (M == MA && N == NB && NA == MB)    /* [MxN] = [MAxNA][MBxNB] */
  1344.     {
  1345.       cblas_dgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
  1346.            alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
  1347.            C->data, INT (C->tda));
  1348.       return GSL_SUCCESS;
  1349.     }
  1350.   else
  1351.     {
  1352.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1353.     }
  1354. }
  1355.  
  1356.  
  1357. int
  1358. gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
  1359.         const gsl_complex_float alpha,
  1360.         const gsl_matrix_complex_float * A,
  1361.         const gsl_matrix_complex_float * B,
  1362.         const gsl_complex_float beta, gsl_matrix_complex_float * C)
  1363. {
  1364.   const size_t M = C->size1;
  1365.   const size_t N = C->size2;
  1366.   const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  1367.   const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  1368.   const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  1369.   const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
  1370.  
  1371.   if (M == MA && N == NB && NA == MB)    /* [MxN] = [MAxNA][MBxNB] */
  1372.     {
  1373.       cblas_cgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
  1374.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1375.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1376.            INT (C->tda));
  1377.       return GSL_SUCCESS;
  1378.     }
  1379.   else
  1380.     {
  1381.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1382.     }
  1383. }
  1384.  
  1385.  
  1386. int
  1387. gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
  1388.         const gsl_complex alpha, const gsl_matrix_complex * A,
  1389.         const gsl_matrix_complex * B, const gsl_complex beta,
  1390.         gsl_matrix_complex * C)
  1391. {
  1392.   const size_t M = C->size1;
  1393.   const size_t N = C->size2;
  1394.   const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  1395.   const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  1396.   const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  1397.   const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
  1398.  
  1399.   if (M == MA && N == NB && NA == MB)    /* [MxN] = [MAxNA][MBxNB] */
  1400.     {
  1401.       cblas_zgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
  1402.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1403.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1404.            INT (C->tda));
  1405.       return GSL_SUCCESS;
  1406.     }
  1407.   else
  1408.     {
  1409.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1410.     }
  1411. }
  1412.  
  1413.  
  1414. /* SYMM */
  1415.  
  1416. int
  1417. gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, float alpha,
  1418.         const gsl_matrix_float * A, const gsl_matrix_float * B,
  1419.         float beta, gsl_matrix_float * C)
  1420. {
  1421.   const size_t M = C->size1;
  1422.   const size_t N = C->size2;
  1423.   const size_t MA = A->size1;
  1424.   const size_t NA = A->size2;
  1425.   const size_t MB = B->size1;
  1426.   const size_t NB = B->size2;
  1427.  
  1428.   if (MA != NA)
  1429.     {
  1430.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1431.     }
  1432.  
  1433.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1434.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1435.     {
  1436.       cblas_ssymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
  1437.            A->data, INT (A->tda), B->data, INT (B->tda), beta,
  1438.            C->data, INT (C->tda));
  1439.       return GSL_SUCCESS;
  1440.     }
  1441.   else
  1442.     {
  1443.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1444.     }
  1445.  
  1446. }
  1447.  
  1448.  
  1449. int
  1450. gsl_blas_dsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, double alpha,
  1451.         const gsl_matrix * A, const gsl_matrix * B, double beta,
  1452.         gsl_matrix * C)
  1453. {
  1454.   const size_t M = C->size1;
  1455.   const size_t N = C->size2;
  1456.   const size_t MA = A->size1;
  1457.   const size_t NA = A->size2;
  1458.   const size_t MB = B->size1;
  1459.   const size_t NB = B->size2;
  1460.  
  1461.   if (MA != NA)
  1462.     {
  1463.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1464.     }
  1465.  
  1466.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1467.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1468.     {
  1469.       cblas_dsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
  1470.            A->data, INT (A->tda), B->data, INT (B->tda), beta,
  1471.            C->data, INT (C->tda));
  1472.       return GSL_SUCCESS;
  1473.     }
  1474.   else
  1475.     {
  1476.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1477.     }
  1478. }
  1479.  
  1480.  
  1481. int
  1482. gsl_blas_csymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1483.         const gsl_complex_float alpha,
  1484.         const gsl_matrix_complex_float * A,
  1485.         const gsl_matrix_complex_float * B,
  1486.         const gsl_complex_float beta, gsl_matrix_complex_float * C)
  1487. {
  1488.   const size_t M = C->size1;
  1489.   const size_t N = C->size2;
  1490.   const size_t MA = A->size1;
  1491.   const size_t NA = A->size2;
  1492.   const size_t MB = B->size1;
  1493.   const size_t NB = B->size2;
  1494.  
  1495.   if (MA != NA)
  1496.     {
  1497.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1498.     }
  1499.  
  1500.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1501.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1502.     {
  1503.       cblas_csymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
  1504.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1505.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1506.            INT (C->tda));
  1507.       return GSL_SUCCESS;
  1508.     }
  1509.   else
  1510.     {
  1511.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1512.     }
  1513. }
  1514.  
  1515. int
  1516. gsl_blas_zsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1517.         const gsl_complex alpha, const gsl_matrix_complex * A,
  1518.         const gsl_matrix_complex * B, const gsl_complex beta,
  1519.         gsl_matrix_complex * C)
  1520. {
  1521.   const size_t M = C->size1;
  1522.   const size_t N = C->size2;
  1523.   const size_t MA = A->size1;
  1524.   const size_t NA = A->size2;
  1525.   const size_t MB = B->size1;
  1526.   const size_t NB = B->size2;
  1527.  
  1528.   if (MA != NA)
  1529.     {
  1530.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1531.     }
  1532.  
  1533.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1534.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1535.     {
  1536.       cblas_zsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
  1537.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1538.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1539.            INT (C->tda));
  1540.       return GSL_SUCCESS;
  1541.     }
  1542.   else
  1543.     {
  1544.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1545.     }
  1546. }
  1547.  
  1548.  
  1549. /* HEMM */
  1550.  
  1551. int
  1552. gsl_blas_chemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1553.         const gsl_complex_float alpha,
  1554.         const gsl_matrix_complex_float * A,
  1555.         const gsl_matrix_complex_float * B,
  1556.         const gsl_complex_float beta, gsl_matrix_complex_float * C)
  1557. {
  1558.   const size_t M = C->size1;
  1559.   const size_t N = C->size2;
  1560.   const size_t MA = A->size1;
  1561.   const size_t NA = A->size2;
  1562.   const size_t MB = B->size1;
  1563.   const size_t NB = B->size2;
  1564.  
  1565.   if (MA != NA)
  1566.     {
  1567.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1568.     }
  1569.  
  1570.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1571.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1572.     {
  1573.       cblas_chemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
  1574.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1575.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1576.            INT (C->tda));
  1577.       return GSL_SUCCESS;
  1578.     }
  1579.   else
  1580.     {
  1581.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1582.     }
  1583.  
  1584. }
  1585.  
  1586.  
  1587. int
  1588. gsl_blas_zhemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1589.         const gsl_complex alpha, const gsl_matrix_complex * A,
  1590.         const gsl_matrix_complex * B, const gsl_complex beta,
  1591.         gsl_matrix_complex * C)
  1592. {
  1593.   const size_t M = C->size1;
  1594.   const size_t N = C->size2;
  1595.   const size_t MA = A->size1;
  1596.   const size_t NA = A->size2;
  1597.   const size_t MB = B->size1;
  1598.   const size_t NB = B->size2;
  1599.  
  1600.   if (MA != NA)
  1601.     {
  1602.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1603.     }
  1604.  
  1605.   if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
  1606.       || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
  1607.     {
  1608.       cblas_zhemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
  1609.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1610.            INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
  1611.            INT (C->tda));
  1612.       return GSL_SUCCESS;
  1613.     }
  1614.   else
  1615.     {
  1616.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1617.     }
  1618. }
  1619.  
  1620. /* SYRK */
  1621.  
  1622. int
  1623. gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
  1624.         const gsl_matrix_float * A, float beta, gsl_matrix_float * C)
  1625. {
  1626.   const size_t M = C->size1;
  1627.   const size_t N = C->size2;
  1628.   const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1629.  
  1630.   if (M != N)
  1631.     {
  1632.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1633.     }
  1634.   else if (N != K)
  1635.     {
  1636.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1637.     }
  1638.  
  1639.   cblas_ssyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
  1640.            INT (A->tda), beta, C->data, INT (C->tda));
  1641.   return GSL_SUCCESS;
  1642. }
  1643.  
  1644.  
  1645. int
  1646. gsl_blas_dsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
  1647.         const gsl_matrix * A, double beta, gsl_matrix * C)
  1648. {
  1649.   const size_t M = C->size1;
  1650.   const size_t N = C->size2;
  1651.   const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1652.  
  1653.   if (M != N)
  1654.     {
  1655.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1656.     }
  1657.   else if (N != K)
  1658.     {
  1659.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1660.     }
  1661.  
  1662.   cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
  1663.            INT (A->tda), beta, C->data, INT (C->tda));
  1664.   return GSL_SUCCESS;
  1665.  
  1666. }
  1667.  
  1668.  
  1669. int
  1670. gsl_blas_csyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1671.         const gsl_complex_float alpha,
  1672.         const gsl_matrix_complex_float * A,
  1673.         const gsl_complex_float beta, gsl_matrix_complex_float * C)
  1674. {
  1675.   const size_t M = C->size1;
  1676.   const size_t N = C->size2;
  1677.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1678.  
  1679.   if (M != N)
  1680.     {
  1681.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1682.     }
  1683.   else if (N != MA)
  1684.     {
  1685.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1686.     }
  1687.  
  1688.   cblas_csyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (MA),
  1689.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
  1690.            GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  1691.   return GSL_SUCCESS;
  1692. }
  1693.  
  1694.  
  1695. int
  1696. gsl_blas_zsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1697.         const gsl_complex alpha, const gsl_matrix_complex * A,
  1698.         const gsl_complex beta, gsl_matrix_complex * C)
  1699. {
  1700.   const size_t M = C->size1;
  1701.   const size_t N = C->size2;
  1702.   const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1703.  
  1704.   if (M != N)
  1705.     {
  1706.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1707.     }
  1708.   else if (N != K)
  1709.     {
  1710.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1711.     }
  1712.  
  1713.   cblas_zsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
  1714.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
  1715.            GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  1716.   return GSL_SUCCESS;
  1717. }
  1718.  
  1719. /* HERK */
  1720.  
  1721. int
  1722. gsl_blas_cherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
  1723.         const gsl_matrix_complex_float * A, float beta,
  1724.         gsl_matrix_complex_float * C)
  1725. {
  1726.   const size_t M = C->size1;
  1727.   const size_t N = C->size2;
  1728.   const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1729.  
  1730.   if (M != N)
  1731.     {
  1732.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1733.     }
  1734.   else if (N != K)
  1735.     {
  1736.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1737.     }
  1738.  
  1739.   cblas_cherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
  1740.            INT (A->tda), beta, C->data, INT (C->tda));
  1741.   return GSL_SUCCESS;
  1742. }
  1743.  
  1744.  
  1745. int
  1746. gsl_blas_zherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
  1747.         const gsl_matrix_complex * A, double beta,
  1748.         gsl_matrix_complex * C)
  1749. {
  1750.   const size_t M = C->size1;
  1751.   const size_t N = C->size2;
  1752.   const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1753.  
  1754.   if (M != N)
  1755.     {
  1756.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1757.     }
  1758.   else if (N != K)
  1759.     {
  1760.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1761.     }
  1762.  
  1763.   cblas_zherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
  1764.            INT (A->tda), beta, C->data, INT (C->tda));
  1765.   return GSL_SUCCESS;
  1766. }
  1767.  
  1768. /* SYR2K */
  1769.  
  1770. int
  1771. gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
  1772.          const gsl_matrix_float * A, const gsl_matrix_float * B,
  1773.          float beta, gsl_matrix_float * C)
  1774. {
  1775.   const size_t M = C->size1;
  1776.   const size_t N = C->size2;
  1777.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1778.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1779.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1780.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1781.  
  1782.   if (M != N)
  1783.     {
  1784.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1785.     }
  1786.   else if (N != MA || N != MB || NA != NB)
  1787.     {
  1788.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1789.     }
  1790.  
  1791.   cblas_ssyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
  1792.         INT (A->tda), B->data, INT (B->tda), beta, C->data,
  1793.         INT (C->tda));
  1794.   return GSL_SUCCESS;
  1795. }
  1796.  
  1797.  
  1798. int
  1799. gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
  1800.          const gsl_matrix * A, const gsl_matrix * B, double beta,
  1801.          gsl_matrix * C)
  1802. {
  1803.   const size_t M = C->size1;
  1804.   const size_t N = C->size2;
  1805.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1806.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1807.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1808.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1809.  
  1810.   if (M != N)
  1811.     {
  1812.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1813.     }
  1814.   else if (N != MA || N != MB || NA != NB)
  1815.     {
  1816.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1817.     }
  1818.  
  1819.   cblas_dsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
  1820.         INT (A->tda), B->data, INT (B->tda), beta, C->data,
  1821.         INT (C->tda));
  1822.   return GSL_SUCCESS;
  1823. }
  1824.  
  1825.  
  1826. int
  1827. gsl_blas_csyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1828.          const gsl_complex_float alpha,
  1829.          const gsl_matrix_complex_float * A,
  1830.          const gsl_matrix_complex_float * B,
  1831.          const gsl_complex_float beta, gsl_matrix_complex_float * C)
  1832. {
  1833.   const size_t M = C->size1;
  1834.   const size_t N = C->size2;
  1835.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1836.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1837.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1838.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1839.  
  1840.   if (M != N)
  1841.     {
  1842.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1843.     }
  1844.   else if (N != MA || N != MB || NA != NB)
  1845.     {
  1846.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1847.     }
  1848.  
  1849.   cblas_csyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
  1850.         GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1851.         INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  1852.   return GSL_SUCCESS;
  1853. }
  1854.  
  1855.  
  1856.  
  1857. int
  1858. gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1859.          const gsl_complex alpha, const gsl_matrix_complex * A,
  1860.          const gsl_matrix_complex * B, const gsl_complex beta,
  1861.          gsl_matrix_complex * C)
  1862. {
  1863.   const size_t M = C->size1;
  1864.   const size_t N = C->size2;
  1865.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1866.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1867.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1868.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1869.  
  1870.   if (M != N)
  1871.     {
  1872.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1873.     }
  1874.   else if (N != MA || N != MB || NA != NB)
  1875.     {
  1876.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1877.     }
  1878.  
  1879.   cblas_zsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
  1880.         GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1881.         INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  1882.   return GSL_SUCCESS;
  1883. }
  1884.  
  1885. /* HER2K */
  1886.  
  1887. int
  1888. gsl_blas_cher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1889.          const gsl_complex_float alpha,
  1890.          const gsl_matrix_complex_float * A,
  1891.          const gsl_matrix_complex_float * B, float beta,
  1892.          gsl_matrix_complex_float * C)
  1893. {
  1894.   const size_t M = C->size1;
  1895.   const size_t N = C->size2;
  1896.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1897.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1898.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1899.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1900.  
  1901.   if (M != N)
  1902.     {
  1903.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1904.     }
  1905.   else if (N != MA || N != MB || NA != NB)
  1906.     {
  1907.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1908.     }
  1909.  
  1910.   cblas_cher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
  1911.         GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1912.         INT (B->tda), beta, C->data, INT (C->tda));
  1913.   return GSL_SUCCESS;
  1914.  
  1915. }
  1916.  
  1917.  
  1918. int
  1919. gsl_blas_zher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
  1920.          const gsl_complex alpha, const gsl_matrix_complex * A,
  1921.          const gsl_matrix_complex * B, double beta,
  1922.          gsl_matrix_complex * C)
  1923. {
  1924.   const size_t M = C->size1;
  1925.   const size_t N = C->size2;
  1926.   const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  1927.   const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  1928.   const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  1929.   const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
  1930.  
  1931.   if (M != N)
  1932.     {
  1933.       GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
  1934.     }
  1935.   else if (N != MA || N != MB || NA != NB)
  1936.     {
  1937.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1938.     }
  1939.  
  1940.   cblas_zher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
  1941.         GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  1942.         INT (B->tda), beta, C->data, INT (C->tda));
  1943.   return GSL_SUCCESS;
  1944.  
  1945. }
  1946.  
  1947. /* TRMM */
  1948.  
  1949. int
  1950. gsl_blas_strmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1951.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
  1952.         const gsl_matrix_float * A, gsl_matrix_float * B)
  1953. {
  1954.   const size_t M = B->size1;
  1955.   const size_t N = B->size2;
  1956.   const size_t MA = A->size1;
  1957.   const size_t NA = A->size2;
  1958.  
  1959.   if (MA != NA)
  1960.     {
  1961.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1962.     }
  1963.  
  1964.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  1965.     {
  1966.       cblas_strmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  1967.            alpha, A->data, INT (A->tda), B->data, INT (B->tda));
  1968.       return GSL_SUCCESS;
  1969.     }
  1970.   else
  1971.     {
  1972.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  1973.     }
  1974. }
  1975.  
  1976.  
  1977. int
  1978. gsl_blas_dtrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  1979.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
  1980.         const gsl_matrix * A, gsl_matrix * B)
  1981. {
  1982.   const size_t M = B->size1;
  1983.   const size_t N = B->size2;
  1984.   const size_t MA = A->size1;
  1985.   const size_t NA = A->size2;
  1986.  
  1987.   if (MA != NA)
  1988.     {
  1989.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  1990.     }
  1991.  
  1992.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  1993.     {
  1994.       cblas_dtrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  1995.            alpha, A->data, INT (A->tda), B->data, INT (B->tda));
  1996.       return GSL_SUCCESS;
  1997.     }
  1998.   else
  1999.     {
  2000.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2001.     }
  2002. }
  2003.  
  2004.  
  2005. int
  2006. gsl_blas_ctrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2007.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
  2008.         const gsl_complex_float alpha,
  2009.         const gsl_matrix_complex_float * A,
  2010.         gsl_matrix_complex_float * B)
  2011. {
  2012.   const size_t M = B->size1;
  2013.   const size_t N = B->size2;
  2014.   const size_t MA = A->size1;
  2015.   const size_t NA = A->size2;
  2016.  
  2017.   if (MA != NA)
  2018.     {
  2019.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2020.     }
  2021.  
  2022.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2023.     {
  2024.       cblas_ctrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2025.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  2026.            INT (B->tda));
  2027.       return GSL_SUCCESS;
  2028.     }
  2029.   else
  2030.     {
  2031.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2032.     }
  2033. }
  2034.  
  2035.  
  2036. int
  2037. gsl_blas_ztrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2038.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
  2039.         const gsl_complex alpha, const gsl_matrix_complex * A,
  2040.         gsl_matrix_complex * B)
  2041. {
  2042.   const size_t M = B->size1;
  2043.   const size_t N = B->size2;
  2044.   const size_t MA = A->size1;
  2045.   const size_t NA = A->size2;
  2046.  
  2047.   if (MA != NA)
  2048.     {
  2049.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2050.     }
  2051.  
  2052.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2053.     {
  2054.       cblas_ztrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2055.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  2056.            INT (B->tda));
  2057.       return GSL_SUCCESS;
  2058.     }
  2059.   else
  2060.     {
  2061.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2062.     }
  2063. }
  2064.  
  2065.  
  2066. /* TRSM */
  2067.  
  2068. int
  2069. gsl_blas_strsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2070.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
  2071.         const gsl_matrix_float * A, gsl_matrix_float * B)
  2072. {
  2073.   const size_t M = B->size1;
  2074.   const size_t N = B->size2;
  2075.   const size_t MA = A->size1;
  2076.   const size_t NA = A->size2;
  2077.  
  2078.   if (MA != NA)
  2079.     {
  2080.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2081.     }
  2082.  
  2083.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2084.     {
  2085.       cblas_strsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2086.            alpha, A->data, INT (A->tda), B->data, INT (B->tda));
  2087.       return GSL_SUCCESS;
  2088.     }
  2089.   else
  2090.     {
  2091.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2092.     }
  2093. }
  2094.  
  2095.  
  2096. int
  2097. gsl_blas_dtrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2098.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
  2099.         const gsl_matrix * A, gsl_matrix * B)
  2100. {
  2101.   const size_t M = B->size1;
  2102.   const size_t N = B->size2;
  2103.   const size_t MA = A->size1;
  2104.   const size_t NA = A->size2;
  2105.  
  2106.   if (MA != NA)
  2107.     {
  2108.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2109.     }
  2110.  
  2111.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2112.     {
  2113.       cblas_dtrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2114.            alpha, A->data, INT (A->tda), B->data, INT (B->tda));
  2115.       return GSL_SUCCESS;
  2116.     }
  2117.   else
  2118.     {
  2119.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2120.     }
  2121. }
  2122.  
  2123.  
  2124. int
  2125. gsl_blas_ctrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2126.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
  2127.         const gsl_complex_float alpha,
  2128.         const gsl_matrix_complex_float * A,
  2129.         gsl_matrix_complex_float * B)
  2130. {
  2131.   const size_t M = B->size1;
  2132.   const size_t N = B->size2;
  2133.   const size_t MA = A->size1;
  2134.   const size_t NA = A->size2;
  2135.  
  2136.   if (MA != NA)
  2137.     {
  2138.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2139.     }
  2140.  
  2141.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2142.     {
  2143.       cblas_ctrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2144.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  2145.            INT (B->tda));
  2146.       return GSL_SUCCESS;
  2147.     }
  2148.   else
  2149.     {
  2150.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2151.     }
  2152. }
  2153.  
  2154.  
  2155. int
  2156. gsl_blas_ztrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
  2157.         CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
  2158.         const gsl_complex alpha, const gsl_matrix_complex * A,
  2159.         gsl_matrix_complex * B)
  2160. {
  2161.   const size_t M = B->size1;
  2162.   const size_t N = B->size2;
  2163.   const size_t MA = A->size1;
  2164.   const size_t NA = A->size2;
  2165.  
  2166.   if (MA != NA)
  2167.     {
  2168.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  2169.     }
  2170.  
  2171.   if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
  2172.     {
  2173.       cblas_ztrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
  2174.            GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
  2175.            INT (B->tda));
  2176.       return GSL_SUCCESS;
  2177.     }
  2178.   else
  2179.     {
  2180.       GSL_ERROR ("invalid length", GSL_EBADLEN);
  2181.     }
  2182. }
  2183.